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Abstract 

Wavelet analysis and compression tools are reviewed and different ap¬ 
plications to study MHD and plasma turbulence are presented. We in¬ 
troduce the continuous and the orthogonal wavelet transform and detail 
several statistical diagnostics based on the wavelet coefficients. We then 
show how to extract coherent structures out of fully developed turbulent 
flows using wavelet-based denoising. Finally some multiscale numerical 
simulation schemes using wavelets are described. Several examples for 
analyzing, compressing and computing one, two and three dimensional 
turbulent MHD or plasma flows are presented. 
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1 Introduction 

Turbulence is ubiquitous and plays a critical role for the plasma stability and 
confinement properties of fusion devices, e.g.^ in the tokamak edge region. Tur¬ 
bulence is a regime of fluid, gas and plasma flows characterized by its highly 
nonlinear dynamics [5]. It exhibits a chaotic, z.e., unpredictible behavior and ro¬ 
tational motion all along a wide range of dynamically active scales. In contrast 
to classical dynamical systems, which are low dimensional and conservative, a 
turbulent flow is a dissipative dynamical system, whose behavior is governed by 
a very large, even maybe infinite, number of degrees of freedom. Each field, e.g., 
velocity, vorticity, magnetic field or current density, strongly fluctuates around 
a mean value and one observes that these fluctuations tend to self-organize into 
so-called coherent structures, ie., vortex tubes in hydrodynamics and vorticty 
sheets and current sheets in magnehydrodynamics (MHD). The presence of co¬ 
herent structures results in the strong spatial and temporal flow intermittency, 
which is a key feature of turbulence. Intermittency is understood here that the 
fluctuations become stronger for decreasing scale and are hence more localized. 
The appropriate tool to study intermittency is the wavelet representation due 
to its intrinsic multiscale nature. Indeed, it yields a sparse multiscale represen¬ 
tation of intermittent fields since wavelets are well localized functions in both 
physical and Eourier space. 

The classical theory of homogeneous turbulence [4] assumes that turbu¬ 
lent flows are statistically stationary and homogeneous. This allows to use the 
Eourier space representation to analyze it (e.^., the energy spectrum is the mod¬ 
ulus of the Eourier transform of the velocity auto-correlation), to model it (e.^.. 
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using Large Eddy Simulation) and to compute it {e.g.^ using spectral methods). 
Hence, since the Fourier representation spreads information among the phases 
of all Fourier coefficients, the structural information (z.e., locality in time and in 
space) is lost when one considers only the modulus of the Fourier coefficients, as 
usually done. This is a major drawback of the classical theory of turbulence and 
the reason why we proposed in 1988 m to replace the Fourier representation by 
the wavelet representation to define new analysis and computational tools able 
to preserve time and space local information. If the Fourier representation is 
well suited to study linear dynamical systems (whose behaviour either persists 
at the initial scale or spreads over larger ones), it is not the case for nonlinear dy¬ 
namical systems for which the superposition principle no more holds (z.e., they 
cannot be decomposed into a sum of independent subsystems to be separately 
studied). Moreover, the evolution of nonlinear dynamical systems develop over 
a wide range of scales, since energy is spread from the initially excited scale 
towards smaller and smaller ones (the so-called energy cascade) until finite-time 
singularities develop (e.^., shocks), unless some dissipative mechanisms damp 
energy and thus avoid its ultra-violet divergence. The art of predicting the 
evolution of nonlinear dynamical systems consists of disentangling their active 
components from their passive components, the former being deterministically 
computed while the latter being discarded or statistically modeled. One thus 
performs a distillation process to only retain the components essential to predict 
the nonlinear behaviour. The wavelet representation is particularly appropriate 
for this since it allows to track the evolution in both space and scale and to only 
retain the degrees of freedom in charge of the nonlinear dynamics. Turbulent 
flows are archetypes of nonlinear dynamical systems and therefore good candi¬ 
dates to be analyzed, modelled and computed using the wavelet representation. 

If we now focus on plasma turbulence we are uneasy about the fact that we 
have two different descriptions, depending on which side of the Fourier transform 
we look from. 

• On the one hand, we have a theory [4] that assumes a nonlinear cascade 
in Fourier space for a range of scales, the so-called ‘inertial range’, where 
the flow kinetic energy is statistically (z.e., for ensemble or time or space 
averages) transferred towards smaller scales until reaching Kolmogorov’s 
scale where molecular dissipation transforms kinetic energy into heat. Un¬ 
der these hypotheses, the theory predicts a power-law behaviour for the 
modulus of the energy spectrum in the inertial range. 

• On the other hand, if we study the flow in physical space we do not have 
yet a predictive theory but only empirical observations (from laboratory 
and numerical experiments) showing the emergence and persistence of 
coherent structures, e.g., blobs and current sheets that concentrate most 
of the kinetic and magnetic energy, even for very high Reynolds number 
flows. 

The classical methods for modeling turbulent flows, e.^.. Large Eddy Sim¬ 
ulation (LES), suppose a scale separation (z.e., a spectral gap) and neglect the 
small scale motions, although their effect onto the large scale motions is sta¬ 
tistically modelled (supposing their dynamics to be linear or slaved to them). 
Unfortunately for those methods we have strong evidence, from both labora¬ 
tory and numerical experiments, that there is no spectral gap since all scales 
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of the inertial range are coupled and nonlinearly interact. Moreover, one ob¬ 
serves that coherent structures play a major dynamical role and are responsible 
for the transport and mixing properties of turbulent flows. In consequence one 
might ask the following questions: are coherent structures the dynamical build¬ 
ing blocks of turbulent flows and can we extract them? If we succeed to do 
so, would it be possible to represent them with a reduced number of degrees of 
freedom and would those be sufficient to compute the flow nonlinear dynamics? 

The aim of this review is to offer a primer on wavelets for both continu¬ 
ous and orthogonal transforms. We then detail different diagnostics based on 
wavelet coefficients to analyze and to compress turbulent flows by extracting 
coherent structures. Examples for experimental data from the tokamak Tore 
Supra (Cadarache, France) and numerical simulation data of resistive drift-wave 
and MHD turbulence, illustrate the wavelet tools. Wavelet-based density esti¬ 
mation techniques to improve particle-in-cell numerical schemes are presented, 
together with a particle-in-wavelet scheme that we developed for solving the 
Vlasov-Poisson equations directly in wavelet space. Coherent Vorticity and 
Current sheet Simulation (CVCS), that applies wavelet filtering to the resis¬ 
tive non-ideal MHD equations, is proposed as a new model for turbulent MHD 
flows. It allows to reduce the number of degrees of freedom necessary to com¬ 
pute them, while capturing the nonlinear dynamics of the flow. This review is 
based on the work and publications we have performed within the last 15 years, 
in collaboration with the CEA-Cadarache and other teams in France, Japan and 
United States. Almost all material presented here has already been published 
in some of our papers (cited in the references), and parts have been adapted 
for this review. Let us only mention few references of wavelet techniques that 
have been used to analyze and quantify plasma turbulence: e.g.^ transients [14], 
bicoherence [I6lll5lllll[l5], intermittency |8] and anisotropy [2| . An exhaustive 
review is beyond the scope of our paper and we focus here exclusively on our 
contributions. 

The outline of this review is the following: first, in section 2 we present 
wavelet analysis tools, including a short primer on continuous and orthogonal 
wavelets. Statistical tools in wavelet coefficient space are also introduced. Sec¬ 
tion 3 focusses on coherent structure extraction using wavelet-based denoising. 
Wavelet-based simulation schemes are reviewed in section 4 and section 5 draws 
some conclusions. 


2 Wavelet analysis 

2.1 Wavelets: a short primer 

2.1.1 Continuous wavelet transform 

The wavelet transform [25| unfolds any signal (e.^., in time) or any field (e.^., 
in three-dimensional space) into both space (or time) and scale (or time scale), 
and possibly directions (for dimensions higher than one). The building block 
of the wavelet transform is the ‘mother wavelet’, 'ip{x) G L^(M) with x G M, 
that is a well-localized function with fast decay at infinity and at least one 
vanishing moment (z.e., zero mean) or more. It is also smooth enough in order 
its Fourier transform, '^(/c), exhibits fast decay for \k\ tending to infinity. From 
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the mother wavelet one then generates a family of wavelets, translated by 6 G M, 
the position parameter, dilated (or contracted) by a G M+, the scale parameter, 
and normalized in L^-norm (ie., ||' 0 a, 6||2 = 1) to obtain the set 

V’a,6(a;) = ^ t/’ (- —(1) 

\ a J 

The wavelet transform of / G (M) is the inner product of / with the analyzing 
wavelets ^a,b^ and the wavelet coefficients, that measure the fluctuations of / 
at scale a and position 6, are 

f{a,b) = {f,i’a,b)= f f{x)'ip* f,{x)dx (2) 

JR 

with * denoting the complex conjugate. The ffinction / is reconstructed from its 
wavelet coefficients, as the inner product of / with the set of analyzing wavelets 

m = -^[ (3) 

Jr+ Jr cl 

where \ijj{k)\‘^k~^dk is a constant that depends on the wavelet ip. Sim¬ 

ilarly to the Fourier transform, the wavelet transform corresponds to a change of 
basis (from physical space to wavelet space) and, since it is an isometry, it pre¬ 
serves the inner product {{f^g) = (f^g)) (Plancherel’s theorem) and conserves 
energy (Parseval’s identity), therefore 

[\fix)\^dx=J^ [ /■ |/(a,6)|2^ (4) 

Jr a 

Note that the wavelet coefficients of the continuous wavelet transform are re¬ 
dundant and therefore correlated. This could be illustrated by the patterns 
one observes within the continuous wavelet coefficients of a white noise, which 
correspond to the correlation between the dilated and translated wavelets (the 
white noise being decorrelated by construction) and visualizes the ‘reproduc¬ 
ing kernel’ of the continuous wavelet transform. Due to the fact that wavelets 
are well localized in physical space, the behaviour of the signal at infinity does 
not play any role. Therefore both wavelet analysis and wavelet synthesis can 
be performed locally, in contrast to the Fourier transform which is intrinsically 
non local (Fourier modes are spread all over space). One can also construct 
peculiar wavelets on a dyadic grid A = (j, i) (z.e., scale is sampled by octaves 
j and space by positions 2~H) that are orthogonal to each other and are used 
to construct wavelet orthonormal bases. In contrast to the continuous wavelet 
coefficients eq. il) that are redundant and correlated, the orthogonal wavelet 
coefficients are decorrelated and non redundant (z.e., a signal sampled on N 
points is perfectly represented by N orthogonal wavelet coefficients only). As 
for the Fourier transform, there exists a Fast Wavelet Transform (FWT) that is 
even faster than the Fast Fourier Transform (FFT) whose operation count for 
a one dimensional signal sampled on N points is proportional to instead of 
A^log 2 N for the FFT. 

2.1.2 Orthogonal wavelet transform 

A discrete wavelet representation is obtained by sampling dyadically the scale 
a and the position b introducing aj = 2~^ and hji = iaj with i^j G Z. The 
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resulting discrete wavelets 


generate orthogonal bases for peculiar wavelets. Figure shows five discrete 
wavelets for j = 3,..., 7 and their corresponding Fourier transforms, the 
modulus \^ji\. Note that the scale 2~^ is related to the wavenumber kj as 

kj = k^2\ ( 6 ) 

where k^ = k\i^{k)\dk/ \ilj{k)\dk is the centroid wavenumber of the cho¬ 

sen wavelet. In Fig. we observe the duality between physical and spectral 
space, namely small scale wavelets are well localized in physical space and badly 
localized in spectral space, and vice-versa. Denoting the support of a wavelet 
in physical space by Ax and the one in spectral space hy Ak the Fourier un¬ 
certainty principle requires that the product Ax Ak is bounded from below. In 



Figure 1: Wavelet representation. Physical space (left) and spectral space 
(right). Note that Ax Ak > C is due to the Fourier uncertainty principle. 

this case the orthogonal wavelet coefficients of a function / G I/^(M) are given 

by 

fji = (7) 

and the corresponding orthogonal wavelet series reads 

fix) = Jjif’jiix)- (8) 

j,iez 

The integral in the continuous reconstruction forumla, eq. (|^, can thus be 
replaced by a discrete sum. In practical applications the infinite sums of the 
wavelet series have to be truncated in both scale and position. Limiting the 
analysis to the largest accessible scale of the domain 2^ = L the scaling function 
associated to the wavelet has to be introduced and the wavelet series becomes 

fix) = f 0Oi(^) ffi ^ ^ (^) 
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where 0 is the scaling function and / = (/, 0oi) the corresponding scaling co¬ 
efficients. The smallest scale 2“"^ is given by the sampling rate of the function 
/ which determines the number of grid points N = 2^. The finite domain size 
implies that the number of positions becomes also finite and, choosing T = 1, 
we obtain the range i = 0,2-^ — 1 for j = 0,J — 1. Figure illustrates 
for an orthogonal spline wavelet the discrete scale-space representation for three 
different scales {j = 6, 7, 8) and positions. There exists a fast wavelet transform 


Wavelets 
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Figure 2: Space-scale representation of an orthogonal spline wavelet at three 
different scales and positions, z.e., '06,6? '07,32, V^spos- The modulus of the Fourier 
transform of three corresponding wavelets is shown in the insert (top, left). 

algorithm which computes the orthogonal wavelet coefficients in 0{N) oper¬ 
ations, therefore even faster than the fast Fourier transform whose operation 
count is 0{Nlog2 N) [29] . 

As example we show in Fig. [^the orthogonal wavelet coefficients of an aca¬ 
demic function presenting discontinuities. We observe that wavelet coefficients 
at small scales only have significant values in the vicinity of the discontinuities. 
Hence only few coefficients are needed to represent the function after discarding 
the small wavelet coefficients. 

Extension to higher dimensions: The orthogonal wavelet representation 
can be extended to represent functions in higher space dimensions using tensor 
product constructions, see e.^., la na EH- Figure shows two-dimensional 
orthogonal wavelets constructed by tensor products. 

The wavelet transform can also be generalized to treat vector-valued func¬ 
tions (e.^., velocity or magnetic fields) in d space dimensions by decomposing 
each vector component into an orthogonal wavelet series. In the following we 
consider a vector field v = for d = 3 sampled at resolution 

N = 2^^^ with periodic boundary conditions. Its orthogonal wavelet series reads 

J-l 7 2^-1 

v{x) = y] y] y] (lo) 

j = 0 ^1=1 21,22,^3=0 
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Figure 3: Academic example: function with two discontinuities and one in its 
derivative (top), corresponding modulus of orthogonal wavelet coefficients in 
logarithmic scale using periodic spline wavelets of degree five. 
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Figure 4: Two-dimensional orthogonal wavelets. Scaling function (top, left) and 
the three associated directional wavelets in the horizontal (top, right), vertical 
(bottom, left) and diagonal (bottom, right) direction. 
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using 3D orthogonal wavelets The basis functions are constructed by 

tensor products of a set of one-dimensional wavelets and scaling functions [9l [29] 
which have been periodized since the boundary conditions considered here are 
periodic. The scale index j varies from 0 to J —1, the spatial index i = (ii, i 2 , '<^ 3 ) 
has 2^^ values for each scale 2~^ and each angle indexed by /i = 1, • • • ,7. The 
three Cartesian directions x = correspond to /i = 1,2,3, while 

/i = 4, 5, 6, 7 denote the remaining diagonal directions. The wavelet coefficients 
measure the fluctuations of v at scale 2~^ and around position 2~H for each of 
the seven possible directions fi. The contribution of the vector field v at scale 
2~^ and direction /i can be reconstructed by summation of Vj^^^i'ipj^^^i{x) over 
all positions i: 

2^-1 

( 11 ) 

042 ,^ 3=0 

The contribution of v at scale 2~^ is obtained by 

7 

'VjiA = ( 12 ) 

fl=l 

For more details on wavelets, we refer the reader to several review articles, 
e.g., [ia|23l|39] and textbooks, e.g., [1129]. 

2.2 Wavelet-based statistical diagnostics 

The physical representation gives access to both position and direction, the 
latter when the space dimension is larger than one. The spectral representation 
gives access to both wavenumber and direction, when the space dimension is 
larger than one, but the information on position is spread among the phases of all 
Fourier coefficients. The wavelet representation combines the advantages of both 
representations, while also giving access to scale. For instance if we consider a 
three-dimensional vector-valued field, its orthogonal wavelet coefficients of each 
of its three components are indexed by three positions, seven directions and 
one scale. Thus using the wavelet representation new statistical diagnostics 
can be designed by computing moments of coefficients using summation, either 
over position, direction or scale, or any combination of them. Second order 
moments correspond to energy distributions (e.^., the energy spectrum), while 
higher order moments allow to compute skewness and flatness. In the following 
we will present scale dependent moments, scale-dependent directional statistics 
and scale dependent topological statistics. By topological statistics we mean 
the statistics of bilinear quantities, like the scalar product of a vector field and 
its curl, e.^., helicity. 

In the following, we give a summary of statistical diagnostics based on 
orthogonal wavelet analysis, here applied to a generic vector field following 
the lines of [35]. Decomposing a vector field into orthogonal wavelets, scale- 
dependent distributions of turbulent flows can be measured, including indiffer¬ 
ent directions and also of different flow components. For example, the energy 
and its spatial fluctuations can be quantified at different length scales and in 
different directions and hence longitudinal or transverse contributions can be 
determined. In the case of an imposed magnetic field the contributions in the 
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directions perpendicular or parallel to it can be distinguished. To this end, sta¬ 
tistical quantities based on the wavelet representation can be introduced, and 
the scale-dependent anisotropy and the corresponding intermittency of MHD 
turbulence can be examined. Here we define intermittency as a departure from 
Gaussianity, which is reflected by increasing flatness when scale decreases. Sand- 
born m introduced this definition in the context of boundary layer flows and 
for a historical overview on intermittency we refer to [38]. Alternative defi¬ 
nitions of intermittency can be found, e.^., in [24], for example a steepening 
of the energy spectrum proposed by Kolmogorov in 1962 [26]. In [28] [27l [42] 
related techniques to quantify the anisotropy of the flow and its intermittency 
have been proposed. They used structure functions of either tensorial compo¬ 
nents or applied the 803 decomposition, which is based on spherical harmonics. 
Structure functions which correspond to moments of increments can be directly 
linked to wavelet decompositions (see, e.^., [38]). The increments are wavelet 
coefficients using the poor man’s wavelet, z.e., the difference of two delta dis¬ 
tributions, which has only one vanishing moment, its mean value. This implies 
that the exponent of the detectable scaling laws is limited by the order of the 
structure function and the scale selectivity is reduced as the frequency localiza¬ 
tion of the poor man’s wavelet is rather bad. These drawbacks can be overcome 
using higher order wavelets. 

2.2.1 Scale dependent moments 

To study the scale-dependent directional statistics we consider the component 

with i = 1,2,3 of a generic vector field v. First we define the g-th order 
moment of the scale-dependent vector Vj{x) = which is here 

either the vector field at scale 2~^ and direction n, or the vector field at 
scale 2“-^ , 

M,[vf] = {{vfy). (13) 

By construction the mean value satisfies = 0 and hence the moments are 

automatically centered. These scale-dependent moments are related to the g-th 
order structure functions, as shown , e.^., in [38]. In the following we consider 
the second order moment which is a scale-dependent quadratic mean 

intensity of and the fourth order moment M 4 ^[vj^^] which contains the scale- 
dependent spatial fluctuations. Both moments are related via the flatness factor. 

In anisotropic turbulence typically a preferred direction can be defined, e.g., 
for low magnetic Reynolds number turbulence, or rotating turbulence. These 
flows have statistical symmetries, which we suppose here with respect to the x^- 
axis. For the remaining perpendicular components, ^=1,2, the average of these 
two components is taken, Mq[vj-] = {Mq[vj^^]^Mq[vj^^]}/2, and the superscript 
T represents the perpendicular contribution. The parallel contribution Vj ^ is 
denoted by . 

The wavelet energy spectrum for is obtained using M 2 ] and eq. ^ , 

( 14 ) 
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where A/cj = {kj— kj)hi2 [3Q[ [T]. It is thus directly related to the Fourier 
energy spectrum and yields a smoothed version [HEO]. The orthogonality of 
the wavelets with respect to scale and direction guarantees that the total energy 
is obtained by direct summation, E = ■ E[v^p] = ^ E[v^p^. 

The standard deviation of the energy spectrum at a given wavenumber kj 
quantifies the spatial variability 



2Aki 


M4vf] - [M 2 [vf]y 


(15) 


The ratio of the fourth and second order moments defines the scale-dependent 
flatness factor, 



(16) 


which quantifies the flow intermittency at scale 2~E 

The scale-dependent flatness is related to the energy spectrum (14) and the 
standard deviation (15), 




2 

+ 1 . 


(17) 


as shown in [6] . This relation illustrates that the spatial variability of the energy 
spectrum is directly reflected by the scale-dependent flatness. 


2.2.2 Scale-dependent directional statistics 


To quantify scale-dependent spatial flow anisotropy and anisotropic flow inter¬ 
mittency we introduce wavelet-based measures. Both component-wise anisotropy 
and directional anisotropy of the flow are considered in the following. For the 


scale-dependent mean energy, E[vj^^]^ the anisotropy measure can be defined 
similarly to the classical Fourier representation. Analoguously this can be ex¬ 
tended for its spatial fluctuations, Using the relation between the scale- 

dependent flatness with the energy spectrum and its spatial fluctuations, eq. 
(17), various measures of anisotropic flow intermittency can be defined. 


Component-wise anisotropy: The scale-dependent component-wise aniso¬ 
tropy is defined by the ratio of perpendicular to parallel energy, and its fluctu¬ 
ation, at a given scale scale 2“-^ , respectively, 

E[v^] (t[v^] 

(^E{kj) = , c„{kj) = -yr- . ( 18 ) 

E[v]] a[v]] 

The scale-dependent mean energy, CE{kj) is a smoothed version of the Fourier 
counterpart c{k). The component-wise anisotropy of the spatial fluctuations is 
quantified by Ca{kj). These measures are directly related to the component-wise 
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flatness factors of v^p^ z.e., F[vp and Fpp as shown in [35]. Combining eqs. 
( [T7| ) and (18) results in 


^ F[v]]-l 


(19) 


which yields a scale-dependent measure of component-wise anisotropic intermit- 
tency. 


Directional anisotropy: Scale-dependent measures for directional anisotropy 
can be defined using ratios of perpendicular to parallel energy and fluctuations 
in longitudinal or transverse directions, 


d^ikj) — 

EKl] 

E[vir 

daikj) 

_ 

d^{kj) = 

EKs] 

dlikj) 

_ 

E[vj-y' 

(T[vj-y 


( 20 ) 

( 21 ) 


The longitudinal direction is denoted by the index L, ie., L = /i = ^. The 
subscript /i = 3 corresponds to a transverse direction of the perpendicular com¬ 
ponents, while T represents the other transverse direction of the perpendicular 
components, z.e., T = /i = 1 for or T = /i = 2 for 

Three principal directions, z.e., fi = 1,2 and 3, out of the seven possible 
directions have been selected for the directional statistics. 

The measures d^{kj) and are smoothed versions of the Fourier rep¬ 

resentation 2e^^\ks)/{e^^\ki) -h e^‘^\k 2 )} and {e^^\ks) + e^‘^\ks)}/{e^^\k 2 ) + 
respectively, following the interpretation of the directional statistics 
proposed in |6]. Furthermore these quantities can be related to second order 
structure functions defined in physical space, and respectively we have: 


2D^^\rh) {D^^\rh) + D^^\rh)} 

{DW{rh) + D(^){ri2)} {Di^Hrh) + D(^){rh)} ' 


Structure functions are defined as the spatial average of velocity increments, 
D^^\r) = + r) — v^^\x)Y). Here consists of contributions of 

to scales larger than 2~F which are obtained by low pass filtering using the 
three-dimensional scaling function at scale 2~F The unit vector of the Cartesian 
direction is denoted hy 1 1 . 

Combining eq. 0 and eqs. 
sures [35] : 

f,L = f 4ikj) 1" _ 1 

^ I 4 (^i)/ F[vl ^]-1 

aT _ f (fej) 1 _ ~ ^ 

“ \dl{kj)f F[vl^]-V 

which quantify the scale-dependent anisotropic intermittency in the transverse 
and longitudinal directions. They measure intermittency, not only in the plane 


(23) 

(24) 


20)-(21), yields directional anisotropy mea- 
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perpendicular or in the direction parallel to for example a magnetic field Bq, 
but also in the longitudinal or transverse directions. These measures are equal 
to one for isotropic fields, and their departure from the value one indicates the 
degree of flow anisotropy. 

2.2.3 Scale-dependent topological statistics 

Considering the velocity field u and the corresponding vorticity uj = V 'x u the 
kinetic helicity, H{x) = u can be defined. The helicity yields a measure of 
the geometrical statistics of turbulence. Integrating the helicity over space one 
obtains the mean helicity H = {u ■ uj). The scale-dependent helicity Hj was 
introduced in [46] and is defined by 


Hj{x) = Uj-u}j 


(25) 


It preserves Galilean invariance, though the kinetic helicity itself does not. 
The corresponding mean helicity is obtained by summing Hj over scale, H = 
(Hj) due to the orthogonality of the wavelet decomposition. 

The relative helicity 


h{x) 


H 

|m| |a)| 


(26) 


defines the cosine of the angle between the velocity and the vorticity at each 
spatial position. The range of h lies between —1 and +1. The scale dependent 
relative helicity can be defined correspondingly 


hj{x) 


Hj 

\Uj\ \0Jj\ 


(27) 


The Euler equations of hydrodynamics conserve the mean kinetic helicity, 

_ (j 

while in ideal MHD turbulence the mean cross helicity H = {u • b) and the 
mean magnetic helicity H = (a • 6) are conserved quantities. Here a is the 
vector potential of the magnetic field b. The scale dependent versions of the 
relative cross and magnetic helicities have been introduced in [48| and are defined 
respectively by 

hf{x) = 


Hf 




u 


(28) 


with H^{x) = u b and 


TT 

hf{x) = ^ 


M 


\\^j\ 


(29) 


with H^{x) = a ■ b. These quantities define the cosine of the angle between 
the two vector fields. 


2.3 Application to 3D MHD turbulence 

In the following we show applications of the above scale-dependent wavelet- 
based measures to three-dimensional incompressible magnetohydro dynamic tur¬ 
bulence. To study the anisotropy we analyze flows with uniformly imposed 
magnetic field considering the quasistatic approximation at moderate Reynolds 
numbers for different interaction parameters [35] . For the geometrical statistics 
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full MHD turbulence without imposed mean field is analyzed [48]. The flows 
are computed by direct numerical simulation with a Fourier pseudo-spectral 
method at resolution 512^ and for further details we refer the reader to the 
respective publications. The flow structure of the quasistatic MHD turbulence 



Figure 5: QS-3D-MHD: Modulus of vorticity for quasistatic 3D MHD at Rx = 
235, with = 0, (left) and N = 2 (right) computed by DNS, from [35] . 

is illustrated in Fig. Shown are isosurfaces of the modulus of vorticity for 
two different interaction parameters N. The interaction parameter characterizes 
the intensity of the imposed magnetic field Hq (here chosen in the z direction) 

(T l 

relative to the flow nonlinearity. It is defined hy N = , where a is the 

electrical conductivity, L the integral length scale, p the density and u' the rms 
velocity. In the case without imposed magnetic field, z.e., N = 0 the flow is 
equivalent to isotropic hydrodynamic turbulence and entangled vortex turbes 
can be observed in Fig. left. For N = 2 the structures are aligned parallel to 
the z direction, ie., the direction of the imposed magnetic field, and the flow is 
thus strongly anisotropic. 

The wavelet energy spectra (Fig. left) yield information on the kinetic 
energy at scale 2~^ and the spatial fluctuations are quantified by the standard 
deviation spectra (Fig. right). All spectra have been multiplied by to 
enhance their differences at small scale. We observe that the spectra decay 
with increasing normalized wavenumber kjp where 77 is the Kolmogorov length 
scale. Furthermore the wavelet spectra (dotted lines) do agree well with the 
corresponding Fourier spectra (solid lines). For larger values of N the spectra 
E[uj-] decay faster for increasing kjp. The standard deviation spectra of uj- 
also decay more rapidly when N becomes larger. 

The scale-dependent anisotropy measures allow to analyze the anisotropy at 
different scales. The scale-dependent component-wise anisotropy CE{kj) shown 
in Fig. left, quantifies the anisotropy of the wavelet mean energy spectrum. 
As expected we find for A" = 0 that CE{kj) ~ 1 as the flow is isotropic. The 
departure from the value one corresponds to flow anisotropy, ie., for values 
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Figure 6 : QS-3D-MHD: Wavelet mean energy spectra (left) E^(kj) together 
with the Fourier energy spectra (solid lines). Wavelet standard deviation spec¬ 
tra (right) (j^{kj). All quantites are shown for the perpendicular velocity 
components. The inset (left) shows the corresponding forcing Fourier spectra 
k^/^Ef{k). From [35] . 


smaller than one the energy of the parallel component is predominant of that of 
the perpendicular component, an obervation which holds for both cases, N = 1 
and N = 2. Furthermore the anisotropy is persistent at the small scales and 
yields smaller values for A" = 2 . Now we examine the anisotropy in different 
directions. Figure]^ right, shows the flow anisotropy of the mean wavelet 
spectrum in the longitudinal direction. We find that this measure yields values 
larger than one for A = 1 and 2 and values close to one for A = 0. For A 7 ^ 0 
the correlation of the velocity component parallel to the imposed magnetic field 
in its longitudinal direction is supposed to be stronger than the correlation of 
the perpendicular components. We also see that the scale dependence gets weak 
for kjf] > 0 . 1 . 




Figure 7: QS-3D-MHD: Component-wise anisotropy measure CE{kj) (left) and 
directional anisotropy measure in the longitudinal direction d^{kj). From [35] . 

The scale-dependent flatness of the perpendicular velocity E[uj-] and of the 

parallel velocity E[u^j], shown in Fig. left, quantify the intermittency of the 
different flow components. In all cases we find that the flatness does indeed 
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increase for decreasing scale. At small scales, kjf] > 1 we also see that the 
flatness is larger for larger values of N. The inset shows that F[u^j] behaves 
similarily. 

The component-wise anisotropy of the intermittency at each scale can be 
quantified with A^{kj)^ see Fig. M right. Again we find that for = 0 values 
close to one are found, as expected due to the isotropy of the flow. For N = 1 
and 2 the component-wise anisotropic intermittency {kj) has values larger 
than one for kjf] > 0.1, which means that the perpendicular velocity becomes 
more intermittent than the parallel velocity at small scales. For N = 2 this 
becomes even more pronounced. 



kjf] 



kjT] 


Figure 8: QS-3D-MHD: Scale-dependent flatness of the perpendicular velocity 
with in the inset the corresponding flatness for the parallel velocity (left). 
Anisotropic measure of intermittency A{kj) (right). From [35]. 

To illustrate the scale-dependent geometric statistics we consider homoge¬ 
neous magnetohydrodynamic turbulence at unit Prandtl number without mean 
magnetic field. The flow has been computed by direct numerical simulation 
at resolution 512^ with random forcing and for further details we refer to [48]. 
Figure shows the PDFs of the relative scale-dependent cross and magnetic 
helicity, and • Figure (left) exhibits two peaks at = ±1 which 
corresponds to a pronounced scale-dependent dynamic alignment. The peaks 
even become larger for smaller scales and thus the probability of alignement 
(or anti-alignement) of the velocity and the magnetic field increases. Figure 
(right) illustrates that the distribution of the scale-dependent magnetic helicity 
becomes more symmetric at small scales. The inset shows that the total rela¬ 
tive magnetic helicity is strongly skewed with a peak at -hi, which is due to the 
presence of substantical mean magnetic helicity. 

3 Extraction of coherent structures using wavelets 

In this section we illustrate the extraction of coherent structures using an al¬ 
gorithm which is based on wavelet denoising. We first describe it for one¬ 
dimensional scalar-valued signals and illustrate its performance on an academic 
test signal. We then generalize the algorithm to higher dimensions and to vector¬ 
valued fields. Finally, different applications to experimental and numerical data 
are shown: 
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Figure 9: 3D-MHD: Scale-dependent PDFs of the relative helicities. Cross 
helicity (left) and magnetic helicity (right). The insets show the PDFs 
of the corresponding total relative helicities From [48]. 


• a scalar-valued signal varying in time measured by a Langmuir probe in 
the scrape-off layer of the tokamak Tore Supra (Cadarache, France), 

• a two-dimensional academic example of the synthetic emissivity of a radi¬ 
ating toric shell with additive noise, 

• experimental movies obtained by a fast camera implemented in Tore Supra, 

• two-dimensional vorticity fields computed for resistive drift-wave turbu¬ 
lence (Hasegawa-Wakatani model) using a pseudo-spectral method, 

• three-dimensional vorticity and current density fields computed for resis¬ 
tive MHD turbulence (incompressible MHD equations) using a pseudo- 
spectral method. 

3.1 Extraction algorithm 

3.1.1 Principle 

We propose a wavelet-based method to extract coherent structures that emerge 
out of turbulent flows, both in fluids (e.^., vortices, shock waves in compressible 
fluids, ...) and in plasmas (e.^., bursts, blobs, ...). The goal is to study their role 
regarding the transport and mixing properties of flows in the turbulent regime. 

For this, we use the wavelet representation that keeps track of both time 
and scale, instead of the Fourier representation that keeps track of frequency 
only. Since there is not yet an universal definition of the coherent structures 
encountered in turbulent flows, we use an apophatic method (introduced in 
Hinduist theology several thousands years ago) where one does not try to define 
what an entity (e.^., a phenomenon, a noumenon, ...) is but rather what it 
is not. We thus agree on the minimal and hopefully consensual statement : 
’coherent structures are not noise\ and propose to define them as : ’coherent 
structures are what remains after denoising’. 

The mathematical definition of noise states that a signal is a noise if it cannot 
be compressed in any functional basis. As a result the shortest description of a 
noise is itself. Note that in most of the cases the experimental noise generated 
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by a measure device does not fit the definition of mathematical noise since it 
could be compressed in at least one functional basis parasite frequencies 

can be removed in the Fourier basis). 

This new way of defining coherent structures allows to process signals and 
fields, but also their cuts or projections (e.^., a probe located at one point 
provides a one dimensional cut of a four dimensional space-time field). Indeed, 
the algorithms commonly used to extract coherent structures cannot work for 
cuts or projections, because they require a template of the structures to extract 
(one would need to take into account how the probe sees all possible translations 
and distortions of the coherent structures). The strength of our algorithm is 
that it treats fields and projections the same way. 

Since we assume that coherent structures are what remains after denoising, 
we need a model, not for the structures themselves, but for the noise. Applying 
‘Ockham’s Razor principle’ (or the daw of parsimony’), we choose as first guess 
the simplest possible model: we suppose the noise to be additive, Gaussian and 
white (z.e., uncorrelated). We then project the turbulent signal (in ID), or 
turbulent field (in higher dimensions), into wavelet space and retain only the 
coefficients having their modulus larger than a given threshold. As threshold 
value we follow Donoho and Johnstone’s proposition of a threshold value that 
depends on the variance of the Gaussian noise we want to remove and on the 
chosen sampling rate m- Since the noise variance is not known a priori for 
turbulent signals (the noise being produced by their intrinsic nonlinear dynam¬ 
ics), we designed a recursive method [3] to estimate it from the variance of the 
weakest wavelet coefficients, z.e., those whose modulus is below the threshold 
value. After applying our algorithm we obtain two orthogonal fields: the coher¬ 
ent field retaining all coherent structures and the incoherent field corresponding 
to the noise. We then check a posteriori that the latter is indeed noise-like (z.e., 
spread all over physical space), Gaussian and uncorrelated (z.e., also spread all 
over Fourier space), and thus confirm the hypotheses we have a priori chosen 
for the noise. 

3.1.2 Wavelet denoising 

We consider a signal s{t) sampled on A" = 2"^ points that we want to denoise, 
assuming the noise to be additive, Gaussian and white. We first project s{t) onto 
an orthogonal wavelet basis and then filter out some of the wavelet coefficients 
thus obtained, Sij. We retain only the wavelet coefficients whose modulus is 
larger than a threshold value. The main difficulty is to estimate it a priori and 
we encounter two possible cases: 

• If we a priori know the noise’s variance cr^, the optimal threshold value is 
given by Donoho and Johnstone’s formula m 

€ = {2a^\nN)^/^ . (30) 

In 1994 they proved m that such a wavelet thresholding method is op¬ 
timal to denoise signals in presence of additive Gaussian white noise, be¬ 
cause it minimizes the maximal L^-error (between the denoised signal and 
the noise-free signal) for functions whose regularity is inhomogeneous, such 
as bursty or intermittent turbulent signals. 
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• If we do not a priori know the variance of the noise, that is the most 
usual case, one should use the wavelet-based recursive algorithm we pro¬ 
posed in 1999 [lais]. This algorithm first estimates the variance of the 
noise by considering the variance of the noisy signal (Jq and computes the 
corresponding threshold 


eo = {2al\nNy/r (31) 

The algorithm splits the wavelet coefficients into two classes: the weak co¬ 
efficients whose modulus is smaller than the threshold, and the remaining 
strong coefficients. It then computes the variance of the weak coefficients 
an to obtain a better estimation of the variance of the noise (estimated 
from the wavelet coefficients using Parseval’s theorem) 

where = {0 < j < J, i = 0,..., — 1} is the index set of the wavelet 

coefficients. The algorithm then replaces eg by = (2cr^ In that 

yields a better estimate of the threshold. This procedure is iterated until 
it reaches the optimal threshold value, when ^ 

In [3] we proved that this algorithm converges for signals having a suf¬ 
ficiently sparse representation in wavelet space, such as the intermittent 
signals encountered in turbulence. We also showed that the larger the sig¬ 
nal to noise ratio is, the faster the convergence. Hence, if the signal s{t) 
is only a noise it converges in one iteration and retains eg as the optimal 
threshold. 

Using the optimal threshold we then separate the wavelet coefficients Sij into 
two contributions: the coherent coefficients whose modulus is larger than e 
and the remaining incoherent coefficients . Finally, the coherent component 
(t) is reconstructed in physical space using the inverse wavelet transform, 
while the incoherent component is obtained as (t) = s{t) — (t). 

3.1.3 Extraction algorithm for one-dimensional signals 

We detail the iterative extraction algorithm for the one-dimensional case and 
quote it from [3]: 

Initialization 

• given the signal s{t) of duration T, sampled on an equidistant grid U = 
iT/N for i = 0, A/^ — 1, with N = 2"^, 

• set n = 0 and perform a wavelet decomposition, ie., apply the Fast 
Wavelet Transform [29] to s to obtain the wavelet coefficients sa for 

• compute the variance a^ of 5 as a rough estimate of the variance of 
the incoherent signal and compute the corresponding threshold eg = 
(2\aNalf^‘^, where (Tq = F 
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• set the number of coefficients considered as noise to Nj = N, z.e., to the 
total number of wavelet coefficients. 


Main loop 

Repeat 


• set = Nj and count the number of wavelet coefficients smaller than 

which yields a new value for 

• compute the new variance from the wavelet coefficiens smaller than 
Cn, i-e., crl+i = f 


l®i*l - 

I 0 else, 


(33) 


and the new threshold = (2 In 
• set n = n + 1 


until {Ni==Nf^). 

Final step 

• reconstruct the coherent signal from the coefficients s^- using the in¬ 
verse Fast Wavelet Transform, where 




( Sji for \Sji\>€n 

\ 0 else 


(34) 


• finally, compute pointwise the incoherent signal (ti) = s{ti) — (ti) for 


End 

Note that the signal is split into s{t) = (t) + (t) and its energy into 

+ cr|, since the coherent and incoherent components are orthogonal, 
z.e., (s^, s^) = 0. 

We use the Fast Wavelet Transform (FWT) [29] that is computed with 
(2MN) multiplications, M being the length of the discrete filter defining the 
orthogonal wavelet used. Remark: for all applications presented in this pa¬ 
per, we use Coifiet 12 wavelets [9|, unless otherwise stated. As long as the 
filter length M < ^ log 2 N, the FWT is faster than the FFT (Fast Fourier 
Transform) computed with N log 2 N operations. Consequently, the extraction 
algorithm requires (2nMN) operations, n being the number of iterations, that 
is small, typically less than log 2 N. 

This algorithm defines a sequence of estimated thresholds (e^)^^!^ and the 
corresponding sequence of estimated variances (cr^)^^j^- In |3| we proved that 
this sequence converges after a finite number of iterations by applying a fixed 
point type argument to the iteration function 


1/2 


^ s,v(^n-|-l) — 


2 In AT 


N 






(35) 
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The algorithm stops after n iterations, when J-s,N{^n) = ^n+i, since the 
number of samples N is finite. In [3] we have also proved that the convergence 
rate depends on the signal to noise ratio {SNR = 10 logio(cr^/c^|)), since the 
smaller the SNR, z.e., the stronger the noise, the faster the convergence is. 
Moreover, if the algorithm is applied to a Gaussian white noise, it converges in 
one iteration only. If it is applied to a signal without noise, the signal is fully 
preserved. In [3] we have also proven the algorithm’s idempotence, z.e., if it 
is applied several times the noise is eliminated the first time and the coherent 
signal will remain the same if the algorithm is reapplied several times. This 
would be the case for a Gaussian filter which, in contrast, is not idempotent. 


3.1.4 Application to an academic test signal 


To illustrate the performance of the iterative algorithm we consider a one¬ 
dimensional noisy test signal s{t) sampled on N = 2^^ = 8192 points (Fig. 10, 
middle). It is made by adding a Gaussian white noise re(t), of mean zero and 
variance = 25, to a piecewise regular academic signal a{t) presenting sev¬ 
eral discontinuities, in the function or in its derivatives (Fig. top). The 
signal to noise ratio is SNR = 10 log]^o(^a/<^w) ~ 11 dB. After applying the 
extraction algorithm we estimate the noise variance to be 25.6 and we obtain 
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a coherent signal {t) very close to the original academic signal a{t) (Fig. 
bottom). The incoherent part {t) is homogeneous and noise like with flatness 
3.03, which corresponds to quasi-Gaussianity. In Fig. (bottom) we observe 
that the coherent signal retains all discontinuities and peaks present in the aca¬ 
demic signal a(t), which is an advantage with respect to standard denoising 
techniques, e.^., low pass Fourier filtering, which have smoothed them. In the 
vicinity of the discontinuities we observe slight overshoots, which are more local 
than the classical Gibbs phenomena and could for example be removed using 
the translation invariant wavelet transform [29]. 


3.1.5 Extension of the algorithm to higher dimensional scalar and 
vector-valued fields 

The extraction algorithm was described in section 3.1.3 for one-dimensional 
scalar-valued signals s{t) varying in time. First, it can be extended to higher¬ 
dimensional scalar fields s{x) varying in space cc G where d is the space 
dimension. To this end the extraction algorithm only requires that the one¬ 
dimensional wavelets are replaced by their equivalent d-dimensional wavelets 
using tensor product constructions, see, e.^., [11291139]. 

Second, the extraction algorithm can also be extended to vector-valued fields 
V = ...,'z;^^^) where each component = 1, ...,d is a scalar valued field. 

The extraction algorithm is then applied to each component of the vector field. 
For thresholding the wavelet coefficients we consider the vector in eq. (10). 


Assuming statistical isotropy of the noise, the modulus of the wavelet coef¬ 
ficient vector is computed. The coherent contribution is then reconstructed 
from those coefficients whose modulus is larger than the threshold defined as 
e = (2/d In A)^/^ where d is the dimensionality of the vector field, a the vari¬ 
ance of the noise and N the total number of grid points. The iterative algorithm 
in section 3.1.3 can then be applied in a straightforward way. 
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Figure 10: Denoising of a piecewise regular signal using iterative wavelet thresh¬ 
olding. Top: original academic signal a{t). Middle: Noisy signal s{t) with a SNR 
= 11 dB. Bottom: Denoised signal (t) with a SNR = 28 dB. 
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To extract coherent structures out of turbulent flows we consider the vor- 
ticity field, which is decomposed in wavelet space. Applying the extraction 
algorithm then yields two orthogonal components, the coherent and incoher¬ 
ent vorticity fields. Subsequently the corresponding induced velocity fields can 
be reconstructed by applying the Biot-Savart kernel, which is the inverse curl 
operator. For MHD turbulence, we consider in addition the current density 
and we likewise split it into two components, the coherent and incoherent cur¬ 
rent density fields. Using Biot-Savart’s kernel we reconstruct the coherent and 
incoherent magnetic fields. 

Note that the employed wavelet bases do not a priori constitute divergence- 
free bases. Thus the resulting coherent and incoherent vector fields are not 
necessarily divergence free. However, we checked that the departure from in¬ 
compressibility only occurs in the dissipative range and remains negligible m- 
Another solution would be to use directly div-free wavelets, but they are much 
more cumbersome to implement (TO]. 


3.2 Application to ID experimental signals from tokamaks 


In [22] we presented a new method to extract coherent bursts from turbulent 
signals. Ion density plasma fluctuations were measured by a fast reciprocating 
Langmuir probe in the scrape-off layer of the tokamak Tore Supra (Cadarache, 
France), for a schematic view we refer to Fig. 11 The resulting turbulent signal 




Figure 11: Left: Position of the reciprocating Langmuir probe in the scrape-off 
layer of the tokamak Tore Supra in Cadarache. Right: Schematic top view of 
the probe. 

is shown in Fig. (top). To extract the coherent burst the wavelet represen- 
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tation is used which keeps track of both time and scale and thus preserves the 
temporal structure of the analyzed signal, in contrast to the Fourier represen¬ 
tation which scrambles it among the phases of all Fourier coefficients. Apply¬ 
ing the extraction algorithm described in section 3.1.3 the turbulent signal in 
Fig. (top) is decomposed into coherent and incoherent components (Fig. 
bottom). Both signals are orthogonal to each other and their properties can 
thus be studied independently. This procedure disentangles the coherent bursts, 
which contain most of the density variance, are intermittent and correlated with 
non-Gaussian statistics, from the incoherent background fluctuations, which are 
much weaker, non-intermittent, noise-like and almost decorrelated with quasi- 
Gaussian statistics. 
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Figure 12: Signal s{t) of duration 8.192 ms, corresponding to the saturation 
current fluctuations measured at 1 MHz in the scrape-off layer of the tokamak 
Tore Supra (Gadarache, France). Top: total signal 5, bottom left coherent part 
Sc, and bottom right incoherent part sj. From [22] . 


The corresponding PDFs are shown in Fig. 13 which confirm that the in¬ 
coherent part is indeed Gaussian like, while the total and coherent signal have 
similar skewed PDFs with algebraic heavy tails for positive signal values. Di¬ 
agnostics based on the wavelet representation were also introduced in [22] 
which allow to compare the statistical properties of the original signals with 
their coherent and incoherent components. The wavelet spectra in comparison 
with classical Fourier spectra (obtained via modified periodograms) in Fig. 14 


(left) confirm that the total and coherent signals have almost the same scale 
energy distribution with a power law behavior close to —5/3. Furthermore the 
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Figure 13: Probability density function p(s) estimated using histograms with 50 
bins. PDF of the total signal s (green dashed line), of the coherent component 
Sc (red solid line) and of the incoherent component sj (blue dotted line, together 
with a Gaussian fit with variance cr| (black dotted line). From [22] . 

wavelet spectra agree well with the Fourier spectra. The incoherent signal yields 
an energy equipartition for more than two magnitudes, which corresponds to 
decorelation in physical space. To quantify the intermittency we plot in Fig. [14] 
(right) the scale dependent flatness of the different signals which shows that 
the coherent contribution extracted from the total signal has the largest values 
at small scale (z.e., high frequency) and is thus the most intermittent. In [22] 
we conjectured that the coherent bursts are responsible for turbulent transport, 
whereas the remaining incoherent fluctuations only contribute to turbulent dif¬ 
fusion. This is confirmed by the resulting energy flux of the total, coherent and 
incoherent parts given in Fig. Note that cross correlation between coherent 



Figure 14: Left: wavelet spectra E{ujj) (lines with symbols) and modified peri- 
odograms Eiuj) (lines) of the total signal s (green and +) , coherent signal sc 
(red and o) and incoherent signal sj (blue and o). Right: corresponding scale 
dependent flatness E vs frequency ujj. The horizontal dotted line E{ujj) = 3 
corresponds to the flatness of a Gaussian process. From [22]. 

and incoherent contributions of the electric potential and the saturation current 
are not shown. 
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Figure 15: Energy flux: total (green), coherent (red) and incoherent (blue). The 
split is made using complex valued wavelets. 

3.3 Application to 2D experimental movies from toka- 
maks 

3.3.1 Tomographic reconstruction using wavelet-vaguelette decom¬ 
position 

Cameras installed in tokamaks aquire images which are difficult to interpret, 
since the three-dimensional structure of the plasma is mapped onto two space 
dimensions and thus flattened in a non-trivial way. This implies that the re¬ 
ceived flux cannot be directly related to the volumic emissivity of the plasma, 
which is a major limitation of such optical diagnostics. The reason is that the 
photons collected by each pixel on the camera sensor have been emitted along 
a corresponding ray, rather than out of a single point in space. Nevertheless 
the three-dimensional radiation can be related to the two-dimensional image 
using tomographic reconstruction, because the dominant structures in tokamak 
edge turbulence happen to be field-aligned filaments, commonly known as blobs. 
They have a higher density than their surroundings, and their structure varies 
more slowly along magnetic field lines than in their orthogonal directions. 

Mathematically the tomographic reconstruction corresponds to an inverse 
problem which has a formal solution under the assumed symmetry, but is ill- 
posed in the presence of noise. Taking advantage of the slow variation of the 
fluctuations along magnetic field lines in tokamaks, this inverse problem can be 
modelled by a helical Abel transform, which is a Volterra integral operator of 
the first kind. In [33] we proposed a tomographic inversion technique, based on 
a wavelet-vaguelette decomposition and coupled with wavelet denoising to ex¬ 
tract coherent structures, that allows to detect individual blobs on the projected 
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movie and to analyse their behaviour. The wavelet-vaguelette decomposition 
(WVD) was introduced by Tchamitchian [43] and used by Donoho [12] to solve 
inverse problems in the presence of localized structures. Tomographic inversion 
using the wavelet-vaguelette decomposition is as an alternative to SVD (Singu¬ 
lar Value Decomposition). SVD and WVD regularize the problem by damping 
the modes of the inverse transform to prevent amplification of the noise, z.e., 
modes below a given threshold are eliminated. For WVD the nonlinear it¬ 
erative thresholding procedure (see section 3.1.3) is applied to the vaguelette 
coefficients. Here Coifiets with two vanishing moments are used [9]. However, 
in contrast to SVD, WVD takes in addition advantage of the spatial localization 
of coherent structures present in the plasma. 

The technicalities of WVD are described in detail in [33] , in the following 
we only explain the principle. The helical Abel transform related the plasma 
light emissivity S (a scalar-valued field) to the integral of the volume emissivity 
received by the camera I = KS, where K is a compact continuous operator. 
The reconstruction of the plasma light emissivity S from I is an inverse problem 
which becomes very difficult when S is corrupted by noise, since computing K~^ 
is an ill-posed problem which amplifies the noise. The vaguelettes are operator 
adapted wavelets and a biorthogonal set of basis functions is obtained from 
the wavelet bases ^px by computing K^px and K'^~^'ipx^ where denotes 

the adjoint inverse operator [43]. Note that vaguelettes inherit the localization 
features of wavelets but may loose the translation and scale invariance, and thus 
the fast wavelet transform cannot be applied anymore. 


3.3.2 Application to an academic example 


To illustrate the method we first consider an academic test case with an given 
emissivity map S, having a uniform radiating shell at constant value one and 


zero elsewhere. A two-dimensional cut in the poloidal plane is shown in Fig. 16 
left. Applying the helical Abel transform we generate the corresponding syn¬ 
thetic image / = KS (Fig. middle). Then we add a Gaussian white noise 


with standard deviation 0.5, which yields the synthetic noisy image (Fig. 16 
right). 



Figure 16: Denoising WVD academic test case with a uniform radiating shell. 
Left: source emission intensity S in the poloidal plane. Middle: corresponding 
noiseless image I = KS in the image plane. Right: noisy image obtained by 
adding Gaussian white noise with variance 0.5. From [33] . 
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Applying the WVD reconstrution to the synthetic noisy image (Fig. 16 


right) gives a denoised emissivity map, a poloidal cut is shown in Fig. left. 
We observe that the main features are preserved, ie., the constant emissivity 
shell is well recovered, besides some spurious oscillations close to discontinuities. 
The corresponding denoised image Id = KSd (Fig.[^ right) illustrates that the 
noise has been successfully removed. A comparison with the standard SVD 
technique in [33] (not shown here) illustrates the superiority of the wavelet- 
vaguelette technique. 
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Figure 17: Denoising WVD academic test case. WVD inversion results. Left: 
reconstructed poloidal emissivity map Sd- Right: denoised image Id = KSd- 
From [33] . 


3.3.3 Application to fast camera data from tokamaks 


Now we present an application to an experimental movie acquired during the 
Tore Supra discharge TS42967, where the plasma was fully detached and stabi¬ 
lized over several seconds using a feedback control. The movie has been obtained 
using a fast camera recording at 40 kHz. Moreover, the time average of the whole 
movie was subtracted from each frame, which helps us to decrease the effect of 
reflection on the chamber wall. The algorithm is then applied directly to the 
fluctuations in the signal instead of the full signal. The experimental conditions 
can be found in [33]. One frame of the movie is shown in Fig. 18, left and 


used as input for the WVD reconstruction algorithm. The resulting emissivity 
map in the poloidal plane, in Fig. [^ middle, shows the presence of localized 
blobs, which propagate counterclockwise as observed in the movies, not shown 
here. Thus their propagation velocity can be determined. The corresponding 
denoised movie frame Id (Fig.[^ right) is obtained by applying the operator K 
to the inverted emissivity map Sd- We observe that the noise has been removed 
and the local features such as blobs and fronts have been extracted. 
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Figure 18: WVD-inversion of a snapshot from a movie obtained from Tore 
Supra, discharge TS42967. Left: noisy frame used as input for the WVD al¬ 
gorithm. Middle: reconstructed emissivity map obtained as a result of WVD. 
Right: denoised frame obtained by applying the operator K to the reconstructed 
emissivity map. From [33] . 

3.4 Application to 2D simulations of resistive drift-wave 
turbulence 

At the edge of the plasma in tokamaks drift-waves play an important role in 
the dynamics and transport. In [7] we considered a two-dimensional slab geom¬ 
etry and performed direct numerical simulations using a two-field model, the 
Hasegawa-Wakatani system which describes the main features of resistive drift- 
waves. The evolution equations for the plasma density fluctuations and the 
electrostatic potential fluctuations are coupled via the adiabaticity parameter 
which models the intensity of the parallel electron resistivity. A Poisson equa¬ 
tion relates the vorticity with the electrostatic potential. The wavelet-based 
coherent vortex extraction method (see section 3.1.3) is then applied in [7] to 
assess the role of coherent vorticity for radial transport and to identify only the 
active degrees of freedom which are responsible for the transport. 

Visualizations of the vorticity field for two regimes, the quasi-hydro dynamic 
case and the quasi-adiabatic case, corresponding respectively to low and high 
collisionality of the plasma, are given in Fig. In both cases coherent vor¬ 
tices can be observed and a dipolar structure is framed by the white rectangles. 
Applying the CVE algorithm we split the vorticity fields into coherent and in¬ 
coherent contributions. In the quasi-hydrodynamics case we find that 1.3% of 
the wavelet coefficients are sufficient to retain 99.9% of the energy, while in 
the quasi-adiabatic case 1.8% of the modes retain 99.0% of the energy. The 
statistical properties of the total, coherent and incoherent vorticity fields are 
assessed in Fig. by plotting the vorticity PDFs and the Fourier enstrophy 
spectra for the two cases. For the quasi-hydro dynamic vorticity the PDFs of 
the total and the coherent field are slightly skewed and exhibit a non-Gaussian 
distribution, while for the quasi-adiabatic case a symmetric almost Gaussian 
like distribution can be observed. The variances of the incoherent parts are 
strongly reduced in both cases with respect to the total fields and the PDFs 
have a Gaussian-like shape. The enstrophy spectra illustrate that coherent and 
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Figure 19: Snapshots of the vorticity field for the quasi-hydro dynamic case (left) 
and for the quasi-adiabatic case (right). Abscissa and ordinate correspond to 
the radial and poloidal position, respectively. The white rectangles indicate the 
selected dipoles. From [7]. 


incoherent contributions exhibit a multiscale behavior. The spectra of total and 
coherent vorticity agree well all over the inertial range. The spectra of the inco¬ 
herent contributions have a powerlaw behavior close to which corresponds to 
an equipartition of kinetic energy. In [7] it is furthermore shown that the radial 
density flux, ie., more than 98%, is indeed carried by these coherent modes. 
In the quasi-hydrodynamic regime, coherent vortices exhibit depletion of the 
polarization-drift nonlinearity as shown in the scatter plot of vorticity against 
the electrostatic potential in Fig. Moreover vorticity strongly dominates 
over strain, in contrast to the quasiadiabatic regime. Details can be found in 

3.5 Application to 3D simulations of resistive MHD tur¬ 
bulence 

In [17] we proposed a method for extracting coherent vorticity sheets and current 
sheets out of three-dimensional homogeneous magnetohydro dynamic (MHD) 
turbulence. To this end the wavelet-based coherent vortex extraction method 
(see section 3.1.3) has been applied to vorticity and current density fields com¬ 
puted by direct numerical simulation (DNS) of forced incompressible MHD tur¬ 
bulence without mean magnetic field at resolution of 512^. Coherent vorticity 
sheets and current sheets are extracted from the DNS data at a given time in¬ 
stant. A visualization of isosurfaces of vorticity and current density of the total, 
coherent and incoherent fields is shown in Fig. The coherent vorticity and 
current density are found to preserve both the vorticity sheets and the current 
sheets present in the total fields while retaining only a few percent of the degrees 
of freedom. The incoherent vorticity and current density are shown to be struc¬ 
tureless and of mainly dissipative nature. The spectral distributions in Fig. 
of kinetic and magnetic energies of the coherent fields only differ in the dissipa¬ 
tive range, while the corresponding incoherent fields exhibit quasi-equipartition 
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Figure 20: Top: PDFs of the vorticity. Bottom: Fourierspectrumoftheenstro- 
phy versus wavenumber. Left: quasi-hydrodynamic case. Right: quasi-adiabatic 
case. Dashed line: total field, solid line: coherent part, dotted line: incoher¬ 
ent part. Note that the coherent contribution (solid) superposes the total field 
(dashed), which is thus hidden under the solid line in all four figures. The 
straight lines indicating power laws are plotted for reference. From [7]. 
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Figure 21: Scatter plot of vorticity against electrostatic potential for the coher¬ 
ent part (top) and incoherent part (bottom). Left: quasi-hydrodynamic case; 
right: quasi-adiabatic case. The red dots correspond to the total field, the blue 
dots correspond to a selected vortex dipole in Fig. 19 From [7]. 
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Figure 22: Isosurfaces of vorticity (left) and current density (right) of the total 
(top), coherent (middle) and incoherent contributions (bottom). From [47] . 
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of energy, corresponding to a slope. The probability distribution functions 
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Figure 23: Kinetic (a) and magnetic (b) energy spectra of the total, coherent and 
incoherent fields. The wavenumber is normalized with the Iroshnikov-Kraichnan 
scale. From [47] . 

(PDFs) of total and coherent fields, for both vorticity and current density, in 
Fig. [^coincide almost perfectly, while the incoherent vorticity and current den¬ 
sity fields have strongly reduced variances. The energy flux shown in Fig. ^E\ 





Figure 24: PDFs of the ^-th component of the velocity (a), vorticity (b), mag¬ 
netic field (c) and current density (d) for the total, coherent and incoherent 
contributions. From m- 

confirms that the nonlinear dynamics is indeed fully captured by the coherent 
fields only. The scale-dependent flatness of the velocity and the magnetic field 
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Figure 25: Contributions to the energy flux normalized by the energy dissipation 
rate versus the wavenumber, which is normalized with the Iroshnikov-Kraichnan 
scale. From m- 


in Fig. [^illustrate that the total and coherent fields have similar scale depen¬ 
dent high order moments and reflect strong intermittency characterized by the 
strong increase of the flatness for decreasing scale. The flatness values of the 
incoherent contributions, of both the velocity and the magnetic field are are 
much smaller and do not increase significanlty for decreasing scale, z.e., they 
are not intermittent. 






Figure 26: Scale-dependent flatness of velocity (a) and magnetic field (b) ver¬ 
sus the wavenumber, which is normalized with the Iroshnikov-Kraichnan scale. 
From [47] . 


4 Wavelet-based simulation schemes 

In the following two wavelet-based methods for solving kinetic plasma equations 
are presented: an application of nonlinear wavelet denoising to improve the con¬ 
vergence of particle-in-cell schemes (PIC) and a particle-in-wavelet scheme for 
solving the Vlasov-Poisson equation directly in wavelet space. We also present 
the Coherent Vorticity and Current sheet Simulation (CVCS) method which 
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extends the Coherent Vorticity Simulation (CVS) [191 [lOl developed for the 
Navier-Stokes equations to the resistive non-ideal MHD equations. Numerical 
examples illustrate the properties and the efficiency of the different methods. 

4.1 Improving particle-in-cell (PIC) schemes by wavelet 
denoising 

Plasma simulations using particles are characterized by the presence of noise, a 
typical feature of Monte-Carlo type simulations. The number of particles, which 
is restricted by the computational resources, limits the statistical sampling and 
thus the accuracy of the reconstructed particle distribution function. 

The discretization error generically known as particle noise due to its random- 
like character of the method quantifies the difference between the distribution 
function reconstructed from a simulation using Np particles and the exact dis¬ 
tribution function. The weak scaling of the error with the number of particles, 
(X 1/however limits the reduction of particle noise by increasing the num¬ 
ber of computational particles in practical applications. This has motivated the 
development of various noise reduction techniques, see, e.^., m, which is of 
importance in the validation and verification of particle codes. 

In [31] we proposed a wavelet-based method for noise reduction in the recon¬ 
struction of particle distribution functions from particle simulation data, called 
wavelet-based density estimation (WBDE). The method was originally intro¬ 
duced in m in the context of statistics to estimate probability densities given 
a finite number of independent measurements. WBDE, as used in m, is based 
on a truncation of the wavelet representation of the Dirac delta function associ¬ 
ated with each particle. The method yields almost optimal results for functions 
with unknown local smoothness without compromising computational efficiency, 
assuming that the particles coordinates are statistically independent. It can be 
viewed as a natural extension of the finite size particles (ESP) approach, with 
the advantage of estimating more accurately distribution functions that have 
localized sharp features. The proposed method preserves the moments of the 
particle distribution function to a good level of accuracy, has no constraints 
on the dimensionality of the system, does not require an a priori selection of a 
global smoothing scale, and is able to adapt locally to the smoothness of the 
density based on the given discrete particle data. Indeed, the projection space is 
determined from the data itself, which allows for a refined representation around 
sharp features, and could make the method more precise than PIC for a given 
computational cost. Moreover, the computational cost of the denoising stage is 
of the same order as one time step of a ESP simulation. 

The underlying idea of WBDE is to expand the sampled particle distribution 
function, represented by a histogram, into an orthogonal wavelet basis using 
the fast wavelet transform. We define the empirical density associated to the 
particles positions Xn for n = 1,..., Np where Np is the number of particles, 

Np 

= W E (36) 

n=l 

and where 5 is the Dirac measure. We then project p^{x) onto an orthogonal 
wavelet basis retaining only scales j such that L < j < J where the scales L 
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and J denote the largest and smallest retained scales, respectively [13]. The 
remaining wavelet coefficients are then thresholded retaining only those whose 
modulus is larger than the scale-dependent threshold where K is a 

constant which depends on the regularity of the solution m- Finally the de- 
noised particle density is obtained by applying an inverse fast wavelet transform. 
In [31] Daubechies wavelets with 6 vanishing moments were used. 

In [H] we treated three cases in order to test how the efficiency of the denois- 
ing algorithm depends on the level of collisionality of the plasma. A strongly 
collisional, weakly collisional and collisionless regimes were considered. For the 
strongly collisional regime we computed particle data of force-free collisional re¬ 
laxation involving energy and pinch-angle scattering. The collisionless regime is 
studied using PIC-data corresponding to bump-on-tail and two-stream instabil¬ 
ities in the Vlasov-Poisson system. The third case of a weakly collisional regime 
is illustrated here using guiding-center particle data of a magnetically confined 
plasma in toroidal geometry. The data was generated with the code DELTA5D. 
Figure]^ shows contour plots of the histogram (top row) and the reconstructed 
densities using WBDE for increasing number of particles. It can be seen that 
the WBDE results in efficiently denoised densities and that the error has been 
reduced by a factor two with respect to the raw histograms as shown in Eig. 




poloidal angle poloidal angle poloidal angle 
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Eigure 27: Contour plots of estimates of Sf for the collisional guiding center 
transport particle data: histogram method (top row) and WBDE method (bot¬ 
tom row). The left, center and right columns correspond to Np = 32 • 10^ (left), 
Np = 128 • 10^ (middle) and Np = 1024 • 10^ (right), respectively. The plots 
show 17 isolines equally spaced within the interval [0.5, 0.5]. Erom [31]. 


4.2 Particle-in-wavelets scheme (PIW) 

In [32] we proposed a new numerical scheme, called particle-in-wavelets, for the 
Vlasov-Poisson equations describing the evolution of the particle distribution 
function / in collisionless plasma, and assessed its efficiency in the simplest case 
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Figure 28: RMS error estimate for collisional guiding center transport particle 
data according to the histogram, the POD, and the wavelet methods. The 
reference density is computed with Np = 1024 • 10^. From m- 


of one spatial dimension. In non-dimensional form the equations read 


dtf + vd^f + dx(t>dyf = 0 

+ 1 - 27r /" f{x,V,t)dv = 0 

Jr 


(37) 

(38) 


where 0 is electric potential. The particle distribution function / is discretized 
using tracer particles, and the charge distribution is reconstructed using wavelet- 
based density estimation (WBDE), discussed in the previous section. The lat¬ 
ter consists in projecting the Delta distributions corresponding to the parti¬ 
cles onto a finite dimensional linear space spanned by a family of wavelets, 
which is chosen adaptively. A wavelet-Galerkin Poisson solver is used to com¬ 
pute the electric potential once the wavelet coefficients of the electron den¬ 
sity p(x, t) = /(x, u, t)dv have been obtained by WBDE. The properties 

of wavelets are exploited for diagonal preconditioning of the linear system in 
wavelet space, which is solved by an iterative method, here conjugated gradi¬ 
ents. Similar to classical PIC codes the interpolation method is compatible with 
the charge assignment scheme. Once the electric field E(x, t) = — t) has 
been interpolated at the particle positions the characteristic trajectories, defined 
by x'{t) = v{t) and v'{t) = —E{x{t),v{t)^t) can be advanced in time using the 
Verlet integrator. 

To demonstrate the validity of the PIW scheme, numerical computations 
of Landau damping and of the two-stream instability have been performed in 
[32] . The stability and accuracy have been assessed with respect to reference 
computations obtained with a precise semi-Lagrangian scheme m- We showed 
that the precision is improved roughly by a factor three compared to a classical 
PIC scheme, for a given number of particles [32], as illustrated in Eig.j^for the 
two-stream instability. We observe that PIW remains uniformely more precise 


for any number of particles thanks to its adaptive properties (Eig. 29, left). 
The total CPU time measured in seconds scaled for the PIW code inversely 
proportional to the number of particles, while for PIC and L-PIW the scaling 
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Figure 29: Comparisons between PIW and PIC for the two-stream instability 
test case. Relative error of the electric field at t = 30, as a function of the 
number of particles (left) and the corresponding computing time (right). Note 
that L-PIW is a variant of PIW where only linear filtering has been applied. 
From [3^ . 


changes when the number of particles is too low for a given spatial resolution. 
However, note that the actual CPU time may depend on the implementation, 
since the PIC code is written in Fortran, while the PIW code is written in C+-h, 
although the same computer was used for both codes. 

4.3 Coherent Vorticity and Current sheet Simulation (CVCS) 

Direct numerical simulation (DNS) of turbulent flow has a large computational 
cost due to the huge number of degrees of freedom to be taken into account. 

The required spatial resolution thus becomes prohibitive, e.^., scaling as 
for hydrodynamics using Kolmogorov like arguments [36]. The CVS method, 
introduced in [191 ED] , proposes to reduce the computational cost by taking only 
into account the degrees of freedom that are nonlinearly active. To this end, 
the coherent structure extraction method (presented in section 3) is combined 
with a deterministic integration of the Navier-Stokes equations. At each time 
step the CVE is applied to retain only the coherent degrees of freedom, typically 
a few percent of the coefficients. Then, a set of neighbor coefficients in space 
and scale, called ‘safety zone’, is added to account for the advection of coherent 
vortices and the generation of small scales due to their interaction. Afterwards 
the Navier-Stokes equations are advanced in time using this reduced set of a 
degrees of freedom. Subsequently, the CVE is applied to reduce the number 
of degrees of freedom and the procedure is repeated for the next time step. 

A graphical illustration, in wavelet coefficient space, of the degrees of freedom 
retained at a given time step, is given in Eig.|^ This procedure allows to track 
the flow evolution in space and scale selecting a reduced number of degrees 
of freedom in a dynamically adaptive way. With respect to simulations on a 
regular grid, much less grid points are used in CVS. 

In [49] we extended CVS to compute 3D incompressible magnetohydro dy¬ 
namic (MHD) turbulent flow and developed a simulation method called coherent 
vorticity and current sheet simulation (CVCS). The idea is to track the time 
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Figure 30: Illustration of the safety zone in wavelet coefficient space used in 
CVS. The degrees of freedom retained by CVE are drawn in red, the adjacent 
coefficients of the safety zone are drawn in green, while the coefficients in blue 
correspond to the inactive degrees of freedom which are not computed. The 
interface 77 , defined in space and scale, separates the region dominated by non¬ 
linear interaction (red) from the region dominated by linear dissipation (blue). 
The horizontal green line corresponds to the Kolmogorov dissipation scale ( 77 ) 
is defined by the statistical mean (either ensemble or space average). 


40 









evolution of both coherent vorticity and coherent current density, z.e., current 
sheets. Both the vorticity and current density fields are, respectively, decom¬ 
posed at each time step into two orthogonal components, corresponding to the 
coherent and incoherent contribution, using an orthogonal wavelet representa¬ 
tion. Each of the coherent fields is reconstructed from the wavelet coefficients 
whose modulus is larger than a threshold, while their incoherent counterparts 
are obtained from the remaining coefficients. The two threshold values depend 
on the instantaneous kinetic and magnetic enstrophies. The induced coherent 
velocity and magnetic fields are computed from the coherent vorticity and cur¬ 
rent density, respectively, using the Biot-Savart kernel. In order to compute 
the flow evolution, one should retain not only the coherent wavelet coefficients 
but also their neighbors in wavelet space, the safety zone. A flowchart sum- 



Figure 31: Flowchart describing the principle of CVCS. The superscripts n and 
n + 1 denote time steps. FWT and FWT“^ denote the fast wavelet transform 
and its inverse. Operators performed in wavelet coefficient space are framed by 
the dashed rectangle. From [49]. 

marizing the principle of CVCS is shown in Fig. [31] and the adaption strategy 
in orthogonal wavelet coefficient space in Fig. 

In [49] CVCS was performed for 3D forced incompressible homogeneous 
MHD turbulence without mean magnetic field, for a magnetic Prandtl num¬ 
ber equal to unity. The Navier-Stokes equations coupled with the induction 
equation were solved with a pseudospectral method using 256^ grid points and 
integrated in time with a Runge-Kutta scheme. Different adaption strategies 
to select the optimal saftey zone for CVCS have been studied. We tested the 
influence of the safety zone and of the threshold, as defined in section 3.1.3, by 
considering three cases: 

• CVCSO with safety zone but without iterating the threshold eo, 

• CVCSl with safety zone but with iterating the threshold once ei. 
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Figure 32: Adaption strategy in wavelet coefficient space used in CVCS: retained 
wavelet coefficients (blue), added wavelet coefficients to ensure a graded tree 
(red) and added wavelet coefficients corresponding to the safety zone (green). 
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• CVCS2 without safety zone but without iterating the threshold eo, 

details can be found in [49]. The quality of CVCS was then assessed by com¬ 
paring the results with a direct numerical simulation. It is found that CVCS 
with the safety zone well preserves the statistical predictability of the turbulent 
flow with a reduced number of degrees of freedom. CVCS was also compared 
with a Fourier truncated simulation using a spectral cutoff filter, where the 
number of retained Fourier modes is similar to the number of the wavelet coef¬ 
ficients retained by CVCSO. Figure [33] shows the percentage of retained wavelet 
coefficients for CVCS (with three different adaption strategies) in comparison 
to Fourier filtering (FTO) with a fixed cut-off wavenumber. The percentage of 
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Figure 33: Evolution of the percentage C of retained wavelet coefficients for 
CVCS with three different adaption strategies in comparison with Fourier fil¬ 
tering (FTO) with a fixed cut-off wavenumber. From [49] . 


CVCSO - 

cvcsi - 

CVCS2 - 

FTO . 



retained kinetic energy, magnetic energy, kinetic enstrophy and magnetic enstro- 
phy for the three different CVCS strategies in comparison with Fourier filtering 


(FTO) is plotted in Fig. 34 


Probability density functions of vorticity and current density, normalized by 
the corresponding standard deviation, in Fig. [35] show that CVCSO and CVCSI 
capture well the high order statistics of the flow, while in FTO and in CVCS2 
the tails of the PDFs are reduced with respect to the DNS results. The energy 
spectra of kinetic and magnetic energy in Fig. confirm that CVCSO and 
CVCSI reproduce perfectly the DNS results in the inertial range, where all 
nonlinear acticity takes place, and only differs in the dissipative range. 

The results thus show that the wavelet representation is more suitable than 
the Fourier representation, especially concerning the probability density func¬ 
tions of vorticity and current density and that only about 13% of the degrees 
of freedom (CVCSO) compared to DNS are sufficient to represent the nonlinear 
dynamics of the flow. A visualization comparing both the vorticity and current 
density field for DNS and CVCSO is presented in Fig. [37| 


5 Conclusion 

This paper reviewed different wavelet techniques and showed several of their ap¬ 
plications to MHD and plasma turbulence. Continuous and orthogonal wavelet 
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Figure 34: Percentage of retained kinetic energy (a), magnetic energy (b), ki¬ 
netic enstrophy (c) and magnetic enstrophy (d) for the three different CVCS 
strategies in comparison with Fourier filtering (FTO). From [49] . 




Figure 35: PDFs of the ^-th component of vorticity (a) and current density (b) 
normalized by the corresponding standard deviation. From m- 
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Figure 36: Kinetic (a) and magnetic energy spectra (b). The wavenumber is 
normalized with the Iroshnikov-Kraichnan scale. From [49]. 


DNS CVCSO 



Figure 37: Visualization of isosurfaces of modulus of vorticity (top) and modulus 
of current density (bottom) for DNS (left) and CVCSO (right). From [49]. 
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transforms were presented and some wavelet-based statistical tools described, af¬ 
ter selecting those most appropriate to study turbulence, such as scale-dependent 
second and higher order moments, intermittency measure, together with scale- 
dependent directional statistical measures. Examples of applications to three- 
dimensional incompressible MHD turbulence, computed by DNS, illustrated 
how the flow intermittency can be quantified and how its anisotropy and he- 
licity might vary with scale. The wavelet-based coherent structure extraction 
algorithm was detailed and validated for a test signal. Different applications 
to experimental and numerical turbulent plasma data, in one, two and three 
dimensions, were shown. The underlying methodology of a wavelet-based to¬ 
mographic reconstruction algorithm for denoising images and movies obtained 
with fast cameras in tokamaks were explained and results were presented. Ap¬ 
plications to an academic example and to fast camera data from Tore Supra 
proved the efficiency of the algorithm to extract blobs and fronts while denoising 
the data. Wavelet-based simulation schemes developed in the context of kinetic 
plasma equations were also described. Results computed with them showed how 
wavelet denoising accelerates the convergence of classical PIC schemes and how 
a particle-in wavelet (PIW) scheme solves the Vlasov-Poisson equation directly 
and efficiently in wavelet space. Concerning the fluid equations, in particular 
the resistive non-ideal MHD equations, the coherent vorticity and current sheet 
simulation (CVCS) methods were explained and examples illustrated the prop¬ 
erties and insights the wavelet-based approach offers in the context of MHD and 
plasma turbulence. 


Acknowledgements 

ME and KS are grateful to Sadri Benkadda for inviting them to give review 
lectures at the ITER International School 2014 “High performance computing 
in fusion science” held in Aix-en-Provence, Erance. The manuscript is based 
on ME’s lecture there, entitled ‘Wavelet transforms and their applications for 
ITER’, given on August 26th, 2014. The authors are also indebted to Wouter 
Bos, Romain Nguyen van yen, Naoya Okamoto and Katsunori Yoshimatsu with 
whom we published several papers on wavelet applications, from which the 
material of this review has been taken. This work was supported by the Erench 
Research Eederation for Eusion Studies carried out within the framework of the 
European Eusion Development Agreement (EEDA). 


References 

[1] P.S. Addison. The Illustrated Wavelet Transform Handbook: Introduetory 
Theory and Applieations in Seienee, Engineering, Medieine and Finanee. 
Taylor & Erancis, London, 2002. 

[2] O. Alexandrova, C. Lacombe and A. Mangeney. Spectra and anisotropy of 
magnetic fluctuations in the Earth’s magnetosheath: Cluster observations. 
Ann. Geophys., 26, 3585-3596, 2008. 


46 



[3] A. Azzalini, M. Farge and K. Schneider. Nonlinear wavelet thresholding : a 
recursive method to determine the optimal threshold value. AppL Comput. 
Harm. Anal, 18(2), 177-185, 2005. 

[4] G.K. Batchelor. The Theory of Homogeneous Turhulenee. Cambridge Uni¬ 
versity Press, 1982. 

[5] D. Biskamp. Nonlinear magnetohydrodynamies. Cambridge University 
Press, 1997. 

[6] W.J.T. Bos, L. Liechtenstein and K. Schneider. Small scale intermittency 
in anisotropic turbulence, Phys. Rev. F, 76, 046310 (2007). 

[7] W.J.T. Bos, S. Futatani, S. Benkadda, M. Farge and K. Schneider. The 
role of coherent vorticity for turbulent transport in resistive drift-wave tur¬ 
bulence. Phys. Plasmas, 15, 072305, 2008. 

[8] V. Carbone, C. Regnoli, E. Martines and V. Antoni. Intermittency and 
self-similarity in plasma edge fluctuations. Phys. Plasmas, 7(2), 445-447, 
2000. 

[9] I. Daubechies. Ten leetures on wavelets. SIAM, Philiadelphia, 1992. 

[10] E. Deriaz, M. Earge and K. Schneider. Craya decomposition using com¬ 
pactly supported biorthogonal wavelets. Appl. Comput. Harm. Anal, 28, 
267-284, 2010. 

[11] D. Donoho and I. Johnstone. Ideal spatial adaptation via wavelet shrinkage. 
Biometrika, 81, 425-455 (1994). 

[12] D. Donoho. Nonlinear solution of linear inverse problems by wavelet- 
vaguelette decomposition. Appl. Comput. Harm. Anal, 2(2), 101-127, 
1995. 

[13] D.L. Donoho, I.M. Johnstone, C. Keryacharian and D. Picard. Density 
estimation by wavelet thresholding. Ann. Statist., 24(2), 508-539, 1996. 

[14] V. Dose, C. Venus and H. Zohm. Wavelet analysis of fusion plasma tran¬ 
sients. Phys. Plasmas, 4(2): 323-328, 1997. 

[15] T. Dudok De Wit, O. Alexandrova, I. Eurno, L. Sorriso-Valvo and C. Zim- 
bardo. Methods for characterising microphysical processes in plasmas. In 
Mierophysies of Cosmie Plasmas, 589-617, Springer, 2014. 

[16] T. Dudok De Wit and V. V. KrasnoselSkikh. Wavelet bicoherence analysis 
of strong plasma turbulence at the Earths quasiparallel bow shock. Phys. 
Plasmas, 2(11): 4307-4311, 1995. 

[17] M. Earge and C. Rabreau. Transformee en ondelettes pour detecter et anal¬ 
yser les structures coherentes dans les ecoulements turbulents bidimension- 
nels. C. R. Aead. Sei. Paris, 307, serie II, 1479-1486, 1988. 

[18] M. Earge. Wavelet transforms and their applications to turbulence. Ann. 
Rev. of Fluid Meeh., 24, 395-457 (1992). 


47 



[19] M. Farge, K. Schneider and N. Kevlahan. Non-Gaussianity and coherent 
vortex simulation for two-dimensional turbulence using adaptive orthonor¬ 
mal wavelet basis. Phys. Fluids, 11, 2187-2201 (1999). 

[20] M. Farge and K. Schneider. Coherent vortex simulation (CVS), a semi- 
deterministic turbulence model using wavelets. Flow, Turb. Combust., 6, 
393-426, 2001. 

[21] M. Farge, K. Schneider, G. Pellegrino, A.A. Wray and R.S. Rogallo. Co¬ 
herent vortex extraction in three-dimensional homogeneous isotropic turbu¬ 
lence: Comparison between CVS and POD decompositions. Phys. Fluids, 
15 (10), 2886-2896, 2003. 

[22] M. Farge, K. Schneider and P. Devynck. Extraction of coherent events 
in turbulent edge plasma using orthogonal wavelets. Phys. Plasmas, 13, 
042304, 2006. 

[23] M. Farge and K. Schneider. Wavelets: applications to turbulence. En¬ 
cyclopedia of Mathematical Physics (Eds. J.-P. Erangoise, G. Naber and 
T.S. Tsun), Elsevier, 408-419, 2006. 

[24] U. Erisch. Turbulence. Cambridge University Press, Cambridge, England, 
1995. 

[25] A. Grossmann and J. Morlet. Decomposition of Hardy functions into square 
integrable wavelets of constant shape. SIAM J. Math. Anal, 15(4), 723- 
736, 1984. 

[26] A. N. Kolmogorov. A refinement of previous hypotheses concerning the 
local structure of turbulence in a viscous incompressible fluid at high 
Reynolds number. J. Eluid Mech., 13, 82-85, 1962. 

[27] S. Kurien and K. R. Sreenivasan. Anisotropic scaling contributions to high- 
order structure functions in high-Reynolds-number turbulence. Phys. Rev. 
E, 62, 2206, 2000. 

[28] S. Kurien and K.R. Sreenivasan. In New Trends in Turbulence, edited by M. 
Lesieur, A. Yaglom, and E. David. EDP Sciences, Les Ulis, Erance, 2001. 

[29] S. Mallat. A wavelet tour of signal processing. Academic Press (1998). 

[30] C. Meneveau. Analysis of turbulence in the orthonormal wavelet represen¬ 
tation. J. Eluid Mech., 232, 469, 1991. 

[31] R. Nguyen van yen, D. del Castilo-Negrete, K. Schneider, M. Earge and 
G.Y. Chen. Wavelet-based density estimation for noise reduction in plasma 
simulation using particles. J. Comput. Phys., 229(8), 2821-2839, 2010. 

[32] R. Nguyen van yen, E. Sonnendriicker, K. Schneider and M. Earge. Particle- 
in-wavelets scheme for the ID Vlasov-Poisson equations. ESAIM: Proceed¬ 
ings, 32, 134-148, 2011. 

[33] R. Nguyen van yen, N. Eedorczak, E. Brochard, G. Bonhomme, K. Schnei¬ 
der, M. Earge and P. Monier-Garbet. Tomographic reconstruction of toka- 
mak plasma light emission from single using wavelet-vaguelette decompo¬ 
sition. Nucl. Fusion, 52, 013005, 2012. 


48 



[34] N. Okamoto, K. Yoshimatsu, K. Schneider and M. Farge. Directional 
and scale-dependent statistics of quasi-static magnetohydrodynamic tur¬ 
bulence. ESAIM: Proceedings, 32, 95-102, 2011. 

[35] N. Okamoto, K. Yoshimatsu, K. Schneider and M. Farge. Small-scale 
anisotropic intermittency in magnetohydro dynamic turbulence at low mag¬ 
netic Reynolds number. Phys. Rev. E, 89, 033013, 2014. 

[36] S. Pope. Turbulent flows. Cambridge University Press, 2000. 

[37] V. A. Sandborn. Measurements of intermittency of turbulent motion in a 
boundary layer. J. Eluid Mech., 6, 221-240, 1959. 

[38] K. Schneider, M. Farge and N. Kevlahan. Spatial intermittency in two- 
dimensional turbulence: a wavelet approach. Woods Hole Mathematics, 
Perspectives in Mathematics and Physics, Vol. 34 (Eds. N. Tongring and 
R.C. Penner), World Scientific, 302-328, 2004. 

[39] K. Schneider and M. Farge. Wavelets: mathematical theory. Encyclopedia 
of Mathematical Physics (Eds. J.-P. Erangoise, G. Naber and T.S. Tsun), 
Elsevier, 426-437, 2006. 

[40] K. Schneider and O. Vasilyev. Wavelet methods in computational fluid 
dynamics. Annu. Rev. Eluid Mech., 42, 473-503, 2010. 

[41] E. Sonnendriicker, J. Roche, P. Bertrand and A. Ghizzo. The semi- 
Lagrangian method for the numerical resolution of the Vlasov equation. 
J. Comp. Phys., 149, 201220, 1999. 

[42] L. Sorriso-Valvo, V. Carbone, R. Bruno and P. Veltri. Persistence of small- 
scale anisotropy of magnetic turbulence as observed in the solar wind. Eu- 
rophys. Lett., 75(5), 832, 2006. 

[43] P. Tchamitchian. Biorthogonalite et theorie des operateurs. Rev. Math. 
Iberoamerica, 3, 163-189, 1987. 

[44] B.P. Van Milligen, E. Sanchez, T. Estrada, C. Hidalgo, B. Branas, B. Car¬ 
reras and L. Garcia. Wavelet bicoherence: a new turbulence analysis tool. 
Phys. Plasmas, 2(8), 3017-3032, 1995. 

[45] B.P. Van Milligen, C. Hidalgo and E. Sanchez. Nonlinear phenomena and 
intermittency in plasma turbulence. Phys. Rev. Lett., 74(3), 395, 1995. 

[46] K. Yoshimatsu, N. Okamoto, K. Schneider, Y. Kaneda and M. Earge. In¬ 
termittency and scale-dependent statistics in fully developed turbulence. 
Phys. Rev. E, 79, 026303, 2009. 

[47] K. Yoshimatsu, Y. Kondo, K. Schneider, N. Okamoto, H. Hagiwara and 
M. Earge. Wavelet based coherent vorticity sheet and current sheet ex¬ 
traction from three-dimensional homogeneous magnetohydro dynamic tur¬ 
bulence. Phys. Plasmas, 16, 082306, 2009. 

[48] K. Yoshimatsu, K. Schneider, N. Okamoto, Y. Kawahara and M. Earge. 
Intermittency and geometrical statistics of three-dimensional homogeneous 
magnetohydrodynamic turbulence: A wavelet viewpoint. Phys. Plasmas, 
18, 092304, 2011. 


49 



[49] K. Yoshimatsu, N. Okamoto, Y. Kawahara, K. Schneider and M. Farge. Co¬ 
herent vorticity and current density simulation of three-dimensional mag¬ 
net ohydro dynamic turbulence using orthogonal wavelets. Geophysical and 
Astrophysical Fluid Dynamics, 107(1-2), 73-92, 2013. 


50 



